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Abstract 

We propose a method for eliminating the truncation error associated with any subspace diagonal- 
ization calculation. The new method, called stochastic error correction, uses Monte Carlo sampling to 
compute the contribution of the remaining basis vectors not included in the initial diagonalization. The 
method is part of a new approach to computational quantum physics which combines both diagonalization 
and Monte Carlo techniques. 

1 Introduction 

In 0] a new approach was proposed for finding the low-energy eigenstates of very large or infinite-dimensional 
quantum Hamiltonians. This proposal combines both diagonalization and Monte Carlo methods, each being 
used to solve a portion of the problem for which the technique is most efficient. The first part of the proposal 
is to diagonalize the Hamiltonian restricted to a subspace containing the most important basis vectors for 
each low energy eigenstate. This may be accomplished either through variational techniques or an ab 
initio method such as quasi-sparse eigenvector (QSE) diagonalization. The second step is to include the 
contribution of the remaining basis vectors by Monte Carlo sampling. The use of diagonalization allows one 
to consider systems with fermion sign oscillations and extract information about wavefunctions and excited 
states. The use of Monte Carlo provides tools to handle the exponential increase in the number of basis 
states for large volume systems. 

The first half of this proposal was discussed in |jj . An adaptive diagonalization algorithm known as the 
quasi-sparse eigenvector method was introduced to find the most important basis vectors for each low energy 
eigenstate. This technique is especially valuable when little is known about the low energy states. It is 
also the only method available which can handle non-orthogonal bases, non-Hermitian Hamiltonians, infinite 
dimensional systems, and which can find several low energy states with like quantum numbers simultaneously. 
In this paper we discuss the second half of the diagonalization/Monte Carlo scheme. We introduce several 
new Monte Carlo techniques which we call stochastic error correction (SEC). There are two general varieties 
of stochastic error correction, methods based on a series expansion and those which are not. The series 
method starts with an eigenvector of the Hamiltonian restricted to some starting subspace and then includes 
the contribution of the remaining basis states as terms in an ordered expansion. The idea is to form a 
perturbative expansion centered around a good non-perturbative starting point. 

As an example of a non-series method we discuss a technique called the stochastic Lanczos method. This 
method again starts with eigenvectors of a Hamiltonian submatrix. Using these as starting vectors, we define 
Krylov vectors, \j) , H \ j) , H 2 \ j) • • • , similar to standard Lanczos diagonalization. The new ingredient is 
that matrix elements between Krylov vectors, (j'\ H n \ j) , are computed using matrix diffusion Monte Carlo. 
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Since the method does not rely on a series expansion, it has the advantage that the starting vectors need 
not be close to the exact eigenvectors. 

One can generate a large class of stochastic error correction methods based on other non-series algorithms, 
various ways of rcsumming the series expansion, or combinations of the two techniques. In this introductory 
paper we concentrate on describing the basic principles and implementation of the series and non-series 
approaches. We also present three test problems which demonstrate the potential of the new approach for a 
range of different problems. In the first example we determine the low energy spectrum of 4> 2+ i theory using 
QSE diagonalization and first order corrections using the series method. In the second example we find 
the low energy spectrum of compact U(l) in 2 + 1 dimensions using the stochastic Lanczos method. In the 
last example we find the ground state of the 2 + 1 dimensional Hubbard model using QSE diagonalization 
and first order series stochastic error correction. In each case we compare with published results in the 
literature. We conclude with a summary and some general comments on the new computational scheme. 



2 Series method 

Let \i) be the eigenvectors of a Hamiltonian H restricted to some subspace S. Let \Aj) be the remaining 
basis vectors in the full space not contained in S. We can represent H as 
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We have used Dirac's bra-ket notation to represent the terms of the matrix. In cases where the basis is 
non-orthogonal or the Hamiltonian is non-Hcrmitian, the precise meaning of terms such as (A\ \ H\l) is the 
action of the dual vector to \Ai) upon the vector H |1). We have written the diagonal terms for the basis 
vectors \Aj) with an explicit factor E for reasons to be explained shortly. 

Let us assume that |1) is close to some exact eigenvector of H which we denote as | lfuii) - More concretely 
we assume that the components of |lf u ii) outside S are small enough so that we can expand in inverse powers 
of the introduced parameter E. In order to simply the expansion we choose to shift the diagonal entries so 
that Ai = 0. 

The series method of stochastic error correction is based on the E" 1 expansion, 
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A fu u = AifT 1 + \"E~ 2 ■■■ . (3) 

It is convenient to choose the normalization of the eigenvector such that the |1) component remains 1. The 
convergence of the expansion is controlled by the proximity of |1) to |lf u ii). If |1) is not at all close to |lf u ii) 
then it will be necessary to use a non-series method such as the stochastic Lanczos method discussed in the 
next section. 

At first order in E~ x we find 
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At second order the contributions are 
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These contributions can be calculated by straightforward Monte Carlo sampling. All that is required is an 
efficient way of generating random basis vectors \Ak) with known probability rates. Let P (Atrial) denote the 
probability of selecting |A tr iai) ° n a given trial. If for example we are calculating the first order correction 
to the eigenvalue, then we have 
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3 Stochastic Lanczos 

We now consider a method called stochastic Lanczos which does not require the starting vectors to be close to 
exact eigenvectors of H . This is essential if the eigenvectors of H are not quasi-sparse and require extremely 
large numbers of basis states to represent accurately. 

Let V be the full Hilbert space for our system. As in the previous section let S be the subspace over 
which we have diagonalized H exactly. Let Pg be the projection operator for S and let Xj and \j) be the 
eigenvalues and eigenvectors of H restricted to S so that 



PsHPs \j) = Xj \j) . 



(11) 



Let Z be an auxiliary subspace, one which contains S but excludes very high-energy states. Let Pz be the 
projection operator for Z. We will choose Z such that PzHPz is bounded above. Let a be a real constant 
which is greater than the midpoint of the minimum and maximum eigenvalues of PzHPz- As n — > oo 
the operator [Pz{H — a)Pz] n maps any given state in Z to the corresponding lowest-energy eigenvector of 
PzHPz with non-zero overlap. 

The stochastic Lanczos method uses the operators [Pz{H — a)Pz] n to approximate the low-energy eigen- 
values and eigenvectors of PzHPz- The goal is to diagonalize H in a subspace spanned by vectors 



\d,j) = [P z (H-a)P z } d \j) 



(12) 
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for several values of d and j. This requires calculating (d',j'\ d,j) and (d',j'\ H \d,j). If our Hamiltonian 
matrix is Hermitian, both of these terms can be written in the general form 

(j'\[P z (H-a)P z r\j). (13) 

Therefore it suffices to determine the matrix 

A n = P s [P Z (H - a)P z ] n P s . (14) 

For non-orthogonal bases and non-Hermitian Hamiltonians, the only change is that we use vectors 

[P z (H-a)P z ] d \j) (15) 

to generate approximate right eigenvectors of H and vectors in the dual space 

(j\[P z (H-a)P z ] d (16) 

to produce approximate left eigenvectors. Adding and subtracting Ps(H — a)Ps, we can rewrite 

A n = P s [P Z (H - a)P z - P S (H - a)P s + P S (H - a)P s ] n P S - (17) 

A n can now be evaluated recursively as 

A n +i — B n+ \ + £ B m {H — a) J 4„- m , (18) 
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where 

B n = Ps [Pz{H - a)P z - P S (H - a)P s ] n P s . (19) 

The components of B n are computed by matrix diffusion Monte Carlo. One could also directly evaluate 
the components of A n . However the calculation for B n eliminates the need to sample the matrix Ps(H—a)Ps, 
which is already known. Any general matrix product M^M^ ■ ■ ■ M^ 71 ' is a sum of degree n monomials, 

M% (2) -mH = £ Mj»Ml%.-M£l k . (20) 

L J J k . . 

We can interpret (^fj|) as a sum over paths through the set of basis vectors of Z, 

|j>-|il>-"--|»n-l>-|fc>, (21) 

with an associated weight M^M^J 2 • ■ ■ ik . The components of B n are sampled using ensembles of 
random walkers. We refer the interested reader to [^| for a review of methods in diffusion Monte Carlo. 

We end the section with a discussion of the fermion sign problem. The sign problem is a general issue 
for any Monte Carlo calculation. For a system with sign oscillations the evaluation of a Euclidean-time 
Green's function involves sums 



with the property that 



Ei \ x i 



(22) 



exp(-c-V-T), (23) 



where V is the volume, T is the Euclidean time, and c is a positive constant. We will refer to this term 
as the cancellation ratio. The exponential dependence on V and T makes computations difficult even for 
small systems. 
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The sign problem will affect the calculation of B n in the stochastic Lanczos algorithm and terms in the 
series method discussed in the previous section. The effect however is different from the sign problem in 
typical Monte Carlo Green's function calculations. Stochastic error correction is a calculation of eigenvalues 
and eigenvectors rather than a sampling of the partition function or the time evolution of a given initial state. 
Therefore the quantity of interest is not exp(-HT) but the action of H or H n on approximate eigenvectors 
of H. Due to homogeneity in H the explicit volume dependence does not appear in the cancellation ratio. 
Instead we find 

~ exp(-fc ■ n), (24) 

Li Pi I 

where k is a positive constant. The sign problem will return if the starting point of the SEC calculation 
is very poor and it becomes necessary to use n such that k ■ n is large. However in many cases k ■ n can 
be kept small even for large n since the most important part of the Hamiltonian, PgHPs, is diagonalized 
exactly. In short the sign problem is less severe because stochastic error correction uses the result of subspace 
diagonalization as its starting point. 




4 (J) 4 theory in 2 + 1 dimensions 

The first example we consider is 4 theory in 2 + 1 dimensions near the 0^—0 symmetry restoration phase 
transition. We will use QSE diagonalization and the series version of stochastic error correction to probe 
the low energy spectrum of the theory on both sides of the phase transition. 

In H Magruder demonstrated the existence of a phase transition in </> 4 +1 by extending Chang's duality 
argument for 4>i+i- The statement of the main result is as follows. Consider the two Lagrange densities 

C+ = \d v 4>d v cl> - - 1 4 + (25) 

= ±cV0d>+ - + (26) 

The counterterm is defined so that in the £+ system the <f> self-energy graphs vanish at zero-momentum 
up to two-loop order. By shifting the field 

(27) 

we note that the same counterterm (S 2 _ (same mass dependence but /i + replaced by fx_) is also sufficient to 
renormalize By equating £ + and £_ we obtain a duality constraint between the two theories. One 

feature of this constraint is that the g — > oo limit of £_ is mapped to the g — > limit of C + . Therefore 
whose reflection symmetry — > —0 is broken at small g, must eventually reach the symmetric phase for 
sufficiently large coupling.]] 

The C- phase transition was studied using quasi-sparse eigenvector diagonalization with stochastic error 
correction. Quantities such as the critical coupling, critical exponents, and the low lying energy spectrum 
were studied and, where possible, compared with Monte Carlo results. A full discussion methods and results 
are presented in || . We will very briefly summarize some of the results below. 

The two spatial dimensions of our system are taken to be a periodic box of size 2L by 2L. We will use 
the modal field formalism to describe the Hamiltonian for the theory.^ In the following we let the vectors 
ft represent ordered integer pairs (n x ,n y ) such that \n x \ , \n y \ < N max . The parameter -/V max corresponds 
with a momentum cutoff scale of A = N max %/L. The modal field Hamiltonian has the form^J 

1 The g — * oo limit of is mapped to the g — > oo limit of £_ and so there is no analogous argument for a phase transition 
in C+. Numerical calculations indicate that there is no phase transition for 
2 We refer the reader to [M for a short introduction. 

3 Counterterms were calculated using finite volume perturbation theory. 
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In Figure 1 we have plotted the lowest energy eigenstates in the rest frame for -/V max = 10 and L = 3tt 
(in units where fx = 1). This choice of parameters corresponds with 441 different momentum modes and 
a momentum cutoff scale of A = 3.33/x. The states shown in Figure 1 are the three lowest eigenstates in 
the even and odd <j) — > — symmetry sectors, and the energies are measured relative to the ground state 
energy. In our calculation QSE diagonalization was used keeping 500 Fock states, and the stochastic error 
correction was computed to first-order using the series method. Error bars shown include statistical error 
and an estimate of the contribution from higher order corrections. We see clear evidence of a second-order 
phase transition near |f = 0.9.Q We have labelled the energies of the states according to their physical 
interpretation in the symmetric phase. E\ is the energy for the one-particle state, E 2 (E 3 ) is for the 
two(three)-particle threshold, and E' 2 (E' 3 ) is for the first state above the two(three)-particle threshold. At 
finite volume these energies are continuous functions of the coupling g. One feature which was also observed 
in <j>i + i is the crossing of energy levels due to the double degeneracy of states in the broken symmetry 
phase. i?3 is connected to a one-particle state in the broken phase while E' 2 is connected to a two-particle 
state. The levels E' 2 and E3 therefore cross near the critical point. 

Another interesting phenomenon is the appearance of a bound state in the broken symmetry phase. In 
both the odd and even symmetry sectors we can measure the ratio of the two-particle to one-particle energies: 

Table 1 



a 

4! 


E' 2 /E 2 


E' 3 /E 3 


0.2 


2.01(4) 


1.98(4) 


0.3 


2.01(4) 


2.05(4) 


0.4 


1.95(4) 


1.96(4) 


0.5 


1.87(4) 


1.87(4) 


0.6 


1.86(4) 


1.82(4) 



These results are consistent with the binding energies reported in [|j and ||, which indicate a ratio of 1.83(3) 
near the critical point. 



5 Compact U(l) in 2+1 dimensions 

Compact U(l) lattice gauge theory in 2 + 1 dimensions is a simple but phenomenologically interesting gauge 
model. It is asymptotically free and in the usual continuum limit describes massless non-interacting photons. 
On the other hand if the continuum limit is reached by rescaling the mass gap to remain constant, one instead 
finds a confining theory of massive bosons E3] . The Hamiltonian has the form 



H = ~E^- 2a: E cos ^- ( 31 ) 



1 A more complete discussion of the critical coupling as well as critical exponents can be found in ji| . 
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In (HI) M are link gauge fields, Op is the sum of the links circuiting a plaquette, 
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and 



e 4 a 
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is the strong coupling parameter, which tends to infinity as the lattice spacing a goes to 0. We follow the 
notation of [|l0| in which an overall constant factor of ^- multiplying the right-hand side of ( |3l"| ) is suppressed. 

The energy levels we measure are therefore in units of — ■ 

The diagonalization of lattice gauge Hamiltonians is constrained by the requirements of gauge invariance. 
To preserve gauge invariance it is most convenient to use a basis which diagonalizes the electric field part of 
the Hamiltonian 



£ 
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As our next example of stochastic error correction we will address the 4x4 lattice system at x = 1 using this 
electric field basis. In it was noted that this poses a challenge to standard diagonalization techniques. 
Even on the small 4x4 lattice a surprisingly large number of states, about 10 7 ~ 10 8 , are needed to 
accurately describe the low energy spectrum at x = 1. This problem can be circumvented by modifying 
the basis states to incorporate more of the physics of the ground state. For example one can introduce a 
disordered background of magnetic flux as suggested in G3], and that approach is followed in an ongoing 



project 15 . However we would like to directly address the problem described in [jll| and show how the 
stochastic Lanczos method handles the proliferation of large numbers of basis states in the original electric 
field basis. 

We will choose our starting subspace S to include all basis states 



\ni>) 



(35) 



which satisfy 



<8, 



(36) 



and which can be reached from the strong coupling vacuum by at most two transitions via the plaquette 
operators exp(±i6>p). We take the auxiliary space Z to be the subspace spanned by basis vectors 



(37) 



Using matrix diffusion Monte Carlo, we diagonalize the subspace formed by the states 

\d,j) = [P z (H-a)P z ] d \j). 



(38) 



j) are the eigenvectors of the Hamiltonian restricted to the original subspace S. In our calculations we use 

2 



a = i lnax and d = 0, 1, • • • 12 for cutoff values 



24, 28,32. 



(39) 



In Table 2 we show the results for the ground state energy E for the different cutoff values L^ ax and the 
extrapolated value at i„ lax = oo. The errors shown are estimated statistical errors. For comparison we 
show the results of [|o) obtained using Green's function Monte Carlo (GFMC). 



7 



Tabic 2 







^max — 28 


-^max — 32 




GFMC 


E 


-7.394(3) 


-7.430(3) 


-7.438(3) 


-7.442(4) 


-7.4432(5) 



In Table 3 we show the masses for the lightest six particles in the system extrapolated to the limit L^ ra = oo. 
We have labelled the particles according to their spin J and sign under conjugation C : A — > —A We also 
include results from Mj for the lowest antisymmetric and symmetric glueballs. 



Table 3 



J c 


Mass 


GFMC 


|0") 


3.03(2) 


3.01(6) 


|0+) 


4.03(3) 


4.05(8) 


|2") 


6.8(1) 




|2+) 


6.8(1) 




|0+) 


7.0(2) 




|0") 


7.1(2) 





The results we find appear in agreement with |10|. Unlike most Monte Carlo algorithms, the SEC method 
is able to find excited states with the same quantum numbers as lower lying states. This was also evident 
in the <^f +1 example where we could track many different states crossing the phase transition. The reason 
for this advantage goes back to the design of stochastic error correction as a Monte Carlo improvement of 
a diagonalization scheme. For the U(l) example one can reliably find the eigenvalues and eigenvectors for 
the first twenty or so states in the low energy spectrum. 



6 Hubbard Model 

The last example we consider is the two-dimensional Hubbard model defined by the Hamiltonian 



H=-t ^2 ( c l<r c ^ + c ^ c ^) + t/ 5Z( c iT Cl T c a c a)- 

<i,j>\ c-=T,J. * 



(40) 



The summation < i,j > is over nearest neighbor pairs. c\ a {cic,) is the creation(annihilation) operator for a 
spin a electron at site i. t is the hopping parameter, and U controls the on-site Coulomb repulsion. The 
model has attracted considerable attention in recent years due to its possible connection to e?-wave pairing 
and stripe correlations in high-T c cuprate superconductors. In spite of its simple form, the computational 
difficulties associated with finding the ground state of the model are substantial even for small systems. 
Fermion sign problems render Monte Carlo simulations ineffective for U positive and away from half-filling, 
and the collective effect of very large numbers of basis Fock states make most diagonalization approaches 
very difficult. A brief overview of the history and literature pertaining to numerical aspects of the Hubbard 
model can be found in fig . 

In terms of momentum space variables, the Hubbard Hamiltonian on an N x N periodic lattice has the 
form 
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(41) 



As a test of our methods, we use QSE diagonalization with stochastic error correction to find the ground 
state energy of the 4x4 Hubbard model with 5 electrons per spin. The corresponding Hilbert space has 
about 2 • 10 7 dimensions. For the QSE diagonalization we use momentum Fock states which diagonalize the 
quadratic part of the Hamiltonian. The Hamiltonian is invariant under the symmetry group generated by 
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reflections about the x and y axes, interchanges between x and y, and interchanges between j and f. We 
find it convenient to work with symmetrized Fock states. We will compute stochastic error corrections to 
first order using the series method. 

In Table 4 we present results for the ground state energy. We encountered no trouble with the sign 
problem, and in fact one can easily see that each term in the first order series expression (|5|) is negative 
definite. The energies are measured relative to the energy of the Fermi sea at U = 0. The errors reported 
are statistical errors associated with the first order SEC calculation. Where available, we compare with the 
results presented in which we label as Exact, Projector Quantum Monte-Carlo (PQMC), and Stochastic 
Diagonalization (SD). Stochastic diagonalization is a subspace diagonalization technique similar to QSE but 
one which uses a different method for selecting the subspace and is based on a variational principle [ p)| ]. 
Although the precise number of basis states used in the SD calculations is not listed, we infer from numbers 
reported for a modified 4x4 Hubbard system that roughly 10 5 states were used.0 



Tabic 4 



Coupling 


States 


QSE 


QSE+SEC 


Exact 


SD 


PQMC 


U = 2t 


100 
500 
1000 


-.4797 
-.4945 
-.5006 


-.50147(5) 
-.50181(3) 
-.50198(1) 


-.50194 


-.5010 


-.44(5) 


u = U 


100 
500 
1000 


-1.620 
-1.748 
-1.800 


-1.8113(4) 
-1.8242(3) 
-1.8302(1) 


-1.8309 


-1.829 


-1.8(2) 


u = m 


500 
1000 
2000 


-2.558 
-2.651 
-2.685 


-2.7073(4) 
-2.7208(2) 
-2.7231(1) 


-2.7245 


-2.723 


-2.9(3) 



Apparently QSE diagonalization with SEC handles the 4x4 system quite well with relatively few states. 
Much larger systems are being studied using both higher series corrections and stochastic Lanczos techniques 




7 Summary and comments 

In this paper we presented two versions of stochastic error correction, the series method and the stochastic 
Lanczos method. The series method starts with eigenvectors of the Hamiltonian restricted to some optimized 
subspace and includes the contribution of the remaining basis states as an ordered expansion. The stochastic 
Lanczos method starts with eigenvectors of a Hamiltonian submatrix and constructs matrix elements of 
Krylov vectors using matrix diffusion Monte Carlo. This method has the advantage that the starting 
vectors need not be close to exact eigenvectors. 

We presented three different examples which demonstrate the potential of the new approach for strongly 
coupled scalar, gauge, and fermionic theories. In the first example we calculated the low energy spectrum 
of 

02+i using the series method, and in the second example we found the spectrum of compact U(l) in 
2 + 1 dimensions using the stochastic Lanczos method. In both examples we found agreement with results 
from the literature. We also found that unlike typical Monte Carlo results, the SEC method is able to find 
the eigenvalues and eigenvectors for excited states with the same quantum numbers as lower lying states. 
This advantage is due to its design as a Monte Carlo improvement of a diagonalization scheme. In the last 
example we found the ground state of the 2 + 1 dimensional Hubbard model using QSE diagonalization and 
first order series stochastic error correction. In this calculation we encountered no fermion sign problem and 
found that our methods yielded very accurate results with far less effort than existing techniques. We believe 
that the methods we have presented hold considerable potential for studying a wide range of non-perturbative 
quantum systems and answering questions difficult to address using other methods. 

Acknowledgments We are grateful to Eugene Golowich, Daniel Lee, and Nikolay Prokofiev for conversa- 
tions. Financial support provided by the National Science Foundation. 

5 The discrete symmetries of the system were not utilized in their calculations. 
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